update_weights Subroutine

    subroutine update_weights(kohonen_map,current_values,ihit,jhit,khit,&
        maximum_radius,iteration) 
!========================================================================================
!!    Subroutine to update the weights   
        class(self_organizing_map) :: kohonen_map
!! A `self_organizing_map` object
        real(kind=wp),dimension(:,:),intent(inout) :: current_values
!! A real array with the values of the current unit
        integer,intent(inout) :: ihit,jhit,khit,iteration
!! Integer variables with the coordinates of the unit (neuron) to be modified
        real(kind=wp),intent(inout) :: maximum_radius
!! Real variable with the maximum radius of the neighborhood 
        real(kind=wp),dimension(size(current_values,1),size(current_values,2)) :: prototype_values
        real(kind=wp),dimension(size(current_values,1),size(current_values,2)) :: winner_values,term1,term2
        integer :: nx,ny,nz,debug_option,ic,current_pos,ineigh,jneigh,kneigh,idbg
        real(kind=wp) :: time_factor,current_radius,alpha,sigma2,h_neighborhood,real_distance,term3
        real(kind=wp) :: distance_ratio,geometric_distance2,eps,current_distance,lambda
        !type(influence_function) :: influence_func
        real(kind=wp),dimension(size(current_values,1),size(current_values,2)) :: v_vector
        real(kind=wp) :: v_vector_norm,r,Psi
        character(len=NUMCHAR) :: m_estimator
!
        nx=kohonen_map%parameters%number_nodes_nx;
        ny=kohonen_map%parameters%number_nodes_ny;
        nz=kohonen_map%parameters%number_nodes_nz;
        debug_option=kohonen_map%parameters%debug_level;
        idbg=kohonen_map%parameters%idbg;
        lambda=2.0_wp*(1.0_wp/maximum_radius);
        time_factor=1.0_wp-dble(iteration)/&
                 dble(kohonen_map%parameters%number_epochs*kohonen_map%parameters%number_patterns);
        !current_radius = max(maximum_radius*real(1001-iteration)/1000.0 + 0.9999999999,4.0d0);
        current_radius = max(maximum_radius*time_factor,4.0_wp);
        !alpha = max(kohonen_map%parameters%learning_rate*(1.0d0-real(iteration)/1000.0),0.01d0);
        alpha = max(kohonen_map%parameters%learning_rate*time_factor,0.01_wp);
        sigma2=current_radius**2;
        !
        m_estimator=trim(kohonen_map%parameters%m_estimator);  
!
        do ic=1,size(kohonen_map%coordinates,1)
            current_pos=position2index(ihit,jhit,khit,nx,ny);
            current_distance=kohonen_map%cells_distances(current_pos,ic)
            if(current_distance < current_radius) then
                geometric_distance2=current_distance**2;
                call index2position(ic,nx,ny,nz,ineigh,jneigh,kneigh);
                !write(*,*) ic,ineigh,jneigh,kneigh,ihit,jhit,khit
                select case(trim(kohonen_map%parameters%neighborhood_type))
                    case('gaussian')
                        h_neighborhood=alpha*dexp(-0.5_wp*geometric_distance2/sigma2);
                    case('bubble')
                        h_neighborhood=alpha;
                end select
                if(debug_option > 0) then
                    write(idbg,*) ihit,jhit,khit,ineigh,jneigh,kneigh
                endif
                select case(trim(kohonen_map%parameters%som_type))
                    case('normal_som')                      
                        call kohonen_map%grid(ineigh,jneigh,kneigh)%get_prototype(prototype_values);
                        prototype_values=prototype_values+h_neighborhood*(current_values-prototype_values);
                        !v_vector=(current_values-prototype_values);
                        !v_vector_norm=dsqrt(sum(v_vector**2));
                        !r=v_vector_norm/sigma;
                        !Psi=influence_func%calculate(m_estimator,r);
                        !prototype_values=prototype_values+sigma*h_neighborhood*Psi*v_vector/v_vector_norm;
                        call kohonen_map%grid(ineigh,jneigh,kneigh)%set_prototype(prototype_values);
                    case('visom')
                        !write(*,*) trim(kohonen_map%parameters%som_type)
                        call kohonen_map%grid(ineigh,jneigh,kneigh)%get_prototype(prototype_values);
                        call kohonen_map%grid(ihit,jhit,khit)%get_prototype(winner_values);
                        real_distance=sum((winner_values-prototype_values)**2);
                        if( (ineigh == ihit) .and. (jneigh == jhit) .and. (kneigh == khit) ) then                           
                             prototype_values=prototype_values+h_neighborhood*(current_values-prototype_values);
                        else
                             distance_ratio=dsqrt(real_distance)/(dsqrt(geometric_distance2)*lambda);
                             term1=(current_values-winner_values);
                             term2=(winner_values-prototype_values);
                             eps=max(1.0_wp*time_factor,0.0_wp);
                             term3=1.0_wp;!((1.0d0-eps)+eps)
                             prototype_values=prototype_values+h_neighborhood*(term1+term2*&
                                         (distance_ratio-1.0_wp)*term3);
                        endif
                        !write(*,*) iteration,dsqrt(real_distance),dsqrt(geometric_distance2)*lambda,distance_ratio
                        call kohonen_map%grid(ineigh,jneigh,kneigh)%set_prototype(prototype_values); 
                    case('robust_som')
                        call kohonen_map%grid(ineigh,jneigh,kneigh)%get_prototype(prototype_values);
                        prototype_values=prototype_values+h_neighborhood*(current_values-prototype_values);
                        ! v_vector=(current_values-prototype_values);
                        ! v_vector_norm=dsqrt(sum(v_vector**2));
                        ! r=v_vector_norm/sigma;
                        ! Psi=influence_func%calculate(m_estimator,r);
                        ! prototype_values=prototype_values+sigma*h_neighborhood*Psi*v_vector/v_vector_norm;
                        call kohonen_map%grid(ineigh,jneigh,kneigh)%set_prototype(prototype_values);
                end select
            endif
        enddo!ic
!
    end subroutine update_weights